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State-of-the-Art  Method  to  Quantify 
Hydrologic  Model  Uncertainty 

by  Brian  E.  Skahill 


PURPOSE:  The  purpose  of  this  document  is  to  report  upon  the  initial  independent 
implementation  and  testing  of  the  Bayesian  Markov  Chain  Monte  Carlo  (MCMC)  algorithm 
Differential  Evolution  Markov  Chain  (DE-MC). 

INTRODUCTION:  To  support  the  US  Army  Corps  of  Engineers’ (US  ACE)  use  of  risk-based 
analysis  in  flood  damage  reduction  studies,  Skahill  (2012)  identified  the  most  promising  state-of- 
the-art  and  practice-oriented  approaches  to  robustly  quantify  hydrologic  and  hydraulic  (H&H) 
model  uncertainty.  Skahill  (2012)  provides  a  path  forward  for  related  work  activities,  including 
software  development,  preparation  of  practice-oriented  guidance  documentation,  and  research  and 
development  directed  at  improving  uncertainty  analysis  algorithm  efficiency. 

Bayesian  Markov  Chain  Monte  Carlo  (MCMC),  and  in  particular  DiffeRential  Evolution  Adaptive 
Metropolis  (DREAM)  (Vrugt  et  al.  2008a,  2009),  and/or  its  basis,  Differential  Evolution  Markov 
Chain  (DE-MC)  (ter  Braak  2006),  was  selected  by  Skahill  (2012)  as  the  state-of-the-art  method  for 
estimating  model  parameter  and  predictive  uncertainty.  The  intent  of  MCMC  is  to  sample  (upon 
completion  of  the  bum-in  period),  via  stochastic  simulation,  from  the  noted  target  equilibrium  (i.e., 
posterior  probability)  distribution.  The  purpose  of  this  document  is  to  report  upon  the  initial 
independent  implementation  and  testing  of  the  Bayesian  MCMC  algorithm  DE-MC.  This  new  DE- 
MC  implementation  differs  notably  from  other  MCMC  implementations  in  that  additional  sampler 
bum-in  (bum-in  is  the  initial  period  when  the  MCMC  sampler  has  not  yet  converged  to  its  target 
equilibrium  distribution)  assessment  heuristics  were  incorporated  into  the  algorithm.  These 
heuristics  attempt  to  support  a  more  robust  assessment  of  sampler  bum-in,  rather  than  solely 
relying  upon  a  quantitative  sampler  convergence  diagnostic  which  can  frequently  prematurely 
misdiagnose  convergence  to  the  equilibrium  target  distribution.  Two  nontrivial  test  cases  serve  as  a 
means  to  verify  the  independent  DE-MC  implementation  and  demonstrate  related  capabilities.  The 
first  test  case  involves  a  bimodal  mixed  normal  target  distribution  and  the  second  test  case  involves 
an  application  of  the  Sacramento  Soil  Moisture  Accounting  (SAC-SMA)  hydrologic  model.  This 
technical  note  concludes  with  brief  remarks  regarding  additional  planned  USACE-ERDC  Bayesian 
MCMC  research  and  development. 

BAYESIAN  MCMC  DE-MC:  Markov  Chain  Monte  Carlo  (MCMC)  is  a  fonnal  Bayesian 
approach  for  estimating  the  posterior  probability  distribution  of  specified  adjustable  model 
parameters.  The  idea  behind  MCMC  is  that  while  one  wants  to  compute  a  probability  density, 
p(6\y),  where  0  and  y  represent  the  vector  of  adjustable  model  parameters  and  the  observed  data, 
respectively,  there  is  the  understanding  that  such  an  endeavor  may  be  impracticable.  Additionally, 
simply  being  able  to  generate  a  large  random  sample  from  the  probability  density  would  be  equally 
sufficient  as  knowing  its  exact  form.  Hence,  the  problem  then  becomes  one  of  effectively  and 
efficiently  generating  a  large  number  of  random  draws  from  p(0\y).  It  was  discovered  that  an 
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efficient  means  to  this  end  is  to  construct  a  Markov  chain,  a  stochastic  process  of  values  that  unfold 
in  time,  with  the  following  properties:  (1)  the  state  space  (set  of  possible  values)  for  the  Markov 
chain  is  the  same  as  that  for  9;  (2)  the  Markov  chain  is  easy  to  simulate  from;  and  (3)  the  Markov 
chain’s  equilibrium  distribution  is  the  desired  probability  density  p(9\y).  By  constructing  such  a 
Markov  chain,  one  could  then  simply  run  it  to  equilibrium  and  subsequently  sample  from  its 
stationary  distribution.  A  Markov  chain  with  the  above  mentioned  properties  can  be  constructed  by 
choosing  a  symmetric  proposal  distribution  and  employing  the  Metropolis  acceptance  probability 
(Metropolis  et  al.  1953)  to  accept  or  reject  candidate  points. 

MCMC  simulation  is  more  efficient  than  other  Monte  Carlo  methods.  The  ability  to  sample  from 
the  posterior  probability  distribution  for  the  specified  adjustable  model  parameters,  p{9\y), 
provides  one  with  the  capacity  to  robustly  address  questions  associated  with  the  deployed 
modeled  scenarios/alternatives  from  a  probabilistic  perspective.  For  example,  a  question  such  as 
“What  is  the  probability  a  specific  simulated  state  variable  such  as  flow  will  be  exceeded?”, 
which  of  course  has  direct  application  not  only  to  flood  risk  management  and  hydrologic  design, 
but  also  to  environmental  and  water  quality  analysis,  among  possible  others,  can  be  robustly 
answered  via  application  of  MCMC.  Moreover,  modeled  sources  of  uncertainty;  viz.,  input 
forcing,  model  parameter,  and  model  structure,  can  all  be  encapsulated  in  0  in  attempts  to 
completely  quantify  model  uncertainty  (Vrugt  et  al.  2008a,  b). 

C.  ter  Braak  (2006)  introduced  Differential  Evolution  Markov  Chain  (DE-MC),  which  combines 
the  salient  features  of  the  global  optimization  method  Differential  Evolution  (DE)  (Stom  and  Price 
1995,  1997)  with  Bayesian  Markov  Chain  Monte  Carlo.  Multiple  chains  are  run  in  parallel  with 
DE-MC  and  learn  from  each  other  by  way  of  jump  proposals  that  are  generated  by  taking  the 
difference  of  two  randomly  selected  chains  from  the  current  population.  The  probability  of 
selecting  the  jump  proposal  is  determined  by  using  the  Metropolis  algorithm  (Metropolis  et  al. 
1953).  The  proposal  vector  for  the  simple  but  effective  DE-MC  algorithm  is  given  as  follows: 


xP  =  x/  +  y  (xri  -  xR2)  +  e  (1) 

where  xp,  x„  y,  xR1,  xR,  and  e  represent  the  proposal  vector  of  dimension  d,  the  zth  of  the  N  chains 
which  constitute  the  evolving  population,  a  weighting  factor,  two  unique  vectors  from  the  current 
population,  excluding  the  zth  chain,  randomly  selected  without  replacement,  and  an  error  term 
sampled  randomly  from  a  symmetric  distribution  with  small  variance  compared  with  that  of  the 
target.  N  and  y  constitute  the  two  parameters  associated  with  the  DE-MC  method.  C.  ter  Braak 
(2006)  suggested  N  =  2d  or  3d  for  simple  unimodal  targets  and  N  =  I  (hi  to  20<7  for  more 
complicated  target  distributions.  Assuming  the  target  distribution  is  multivariate  nonnal,  the 
optimal  choice  for  y  is  2.38/2 d,  and  this  value  was  observed  to  work  well  for  the  tests  and 
examples  considered  by  C.  ter  Braak  (2006).  Moreover,  C.  ter  Braak  (2006)  demonstrated  that 
adapting  the  DE-MC  algorithm  such  that  y  equals  one  every  tenth  generation  mitigates  against 
the  potential  of  becoming  trapped  in  a  single  mode  within  a  multimodal  distribution. 

INDEPENDENT  IMPLEMENTATION  AND  TESTING  OF  DE-MC:  The  author  of  this 
technical  note  wrote  an  independent  implementation  of  the  MCMC  method  DE-MC,  as  documented 
herein.  Noteworthy  elements  of  the  independent  implementation  include  the  following: 

1 .  Treatment  of  the  initialization  of  the  population. 
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2.  Treatment  of  jump  proposal  dimensions  that  are  outside  of  their  bounds. 

3.  Use  of  the  Gelinan  and  Rubin  GR  diagnostic  (Gelinan  and  Rubin  1992),  R,  to  assess 
convergence  to  the  stationary  target  distribution. 

Initialization  of  Population.  The  population  to  be  evolved  can  be  initialized  either  via 
uniform  random  sampling  (URS)  or  Latin  hypercube  sampling  (LHS).  Based  on  limited 
experimentation  to  date,  it  is  recommended  that  one  employ  LHS  to  initialize  the  population. 

Treatment  of  Proposal  Dimensions  that  are  out  of  Bounds.  For  DE-MC,  wherein  jump 
proposals  are  defined  by  Equation  1,  it  is  possible  for  one  or  more  dimensions  associated  with  the 
proposal  of  a  given  chain  to  be  out  of  bounds.  In  that  case,  the  specific  proposal  dimension  which 
is  out  of  bounds  can  either  be  (a)  set  equal  to  its  appropriate  boundary  value,  or  alternatively  it  can 
be  (b)  set  equal  to  its  current  value.  Limited  experience  to  date  suggests  recommending  the  use  of 
the  latter  option  (i.e.,  (b))  for  this  implementation  issue. 

Use  of  the  Gelman  and  Rubin  Convergence  Diagnostic.  A  notable  departure  of  the 
independent  DE-MC  implementation  described  herein  from  previous  hydrologic  modeling  studies, 
which  have  employed  Bayesian  MCMC  for  optimization  and  inference,  is  how  the  Gelman  and 
Rubin  quantitative  diagnostic,  R,  is  employed  to  assess  sampler  convergence;  viz.,  how  R  is  used  to 
detennine  if  the  burn-in  period  for  the  MCMC  sampler  is  complete  in  that  the  chains  have  reached 
equilibrium  and  sampling  is  now  from  the  target  distribution.  The  Gelman  and  Rubin  convergence 
diagnostic  effectively  measures  the  within  and  between  variance  of  the  chains,  which  are 
effectively  the  same  when  the  chains  have  sufficiently  mixed.  A  common  specified  value  for  R  to 
diagnose  sampler  convergence  is  1.2,  and  when  the  threshold  value  is  reached,  for  each  element  of 
9,  it  is  often  assumed  that  the  MCMC  sampler  has  converged  to  the  target  distribution  (Vrugt  et  al. 
2003;  Vrugt  et  al.  2008a,  2008b;  Vrugt  et  al.  2009).  The  interested  reader  is  directed  to  Cowles  and 
Carlin  (1996),  and  references  cited  therein,  for  more  detail  and  discussion  pertaining  to  the  Gelman 
and  Rubin  quantitative  convergence  diagnostic.  Unfortunately,  it  is  known  that  diagnostic  has  the 
tendency  to  prematurely  assess  MCMC  sampler  convergence  to  its  stationary  distribution  (Kass  et 
al.  1998).  This  characteristic  for  R  was  observed  during  initial  testing  and  evaluation  of  the 
independent  DE-MC  implementation. 

In  attempts  to  ensure  a  more  robust  assessment  of  DE-MC  sampler  burn-in,  a  hybrid  semi- 
automated  approach  was  implemented,  consistent  with  available  guidance  regarding  practical 
application  of  Bayesian  MCMC.  Effectively,  additional  heuristics  have  been  incorporated  to 
assist  with  diagnosing  sampler  convergence.  In  particular,  if  the  Gelman  and  Rubin  convergence 
diagnostic  indicates  convergence,  in  that  for  each  adjustable  model  parameter  the  value  for  R  is 
less  than  a  pre-specified  threshold  value,  then  diagnostic  information  associated  with  each  chain 
and  for  each  parameter  is  written  to  file  for  manual  review.  In  addition,  the  MCMC  algorithm  is 
paused,  the  user  is  prompted  that  indicates  sampler  convergence  and  that  data  is  available  for 
evaluation  to  support  determination  as  to  whether  or  not  the  burn-in  period  is  complete.  If  the 
user  decides  that  the  sampler  has  not  converged,  then  the  algorithm  proceeds  as  previously,  but 
with  a  reduced  value  for  R ,  the  new  threshold  value  given  by  (l+Rprev)/2.  Alternatively,  if  the 
user  decides  that  the  DE-MC  sampler  has  converged,  then  subsequent  sampling  during  the  post 
bum-in  monitoring  period  proceeds  unless  it  is  determined  that  convergence  was  identified 
prematurely  based  on  a  comparison  of  the  currently  computed  minimum  RMSE  fitness  value 
versus  the  minimum  root-mean-square  error  (RMSE)  fitness  value  computed  at  what  was 
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previously  thought  to  be  the  completion  of  the  bum-in  phase  for  the  sampler.  In  this  case,  burn-in 
sampling  recommences  with  a  new  value  for  R  as  mentioned  previously.  The  threshold 
comparison  value  of  the  noted  minimum  RMSEs,  which  suggests  that  burn-in  was  previously 
assessed  prematurely,  is  a  specified  user  input.  This  additional  heuristic  incorporated  into  the 
algorithm  to  automatically  self-correct  during  post  bum-in  monitoring  is  relatively  simplistic  and 
requires  further  evaluation  and  testing.  Testing  to  date,  has  demonstrated  improved  assessment  of 
DE-MC  MCMC  sampler  convergence. 

Verification  and  Testing.  Verification  of  the  independent  DE-MC  implementation  involved 
two  documented  tests:  one,  a  non-trivial  known  target  distribution  (ter  Braak  2006),  and  the  other 
involving  calibration  of  a  Sacramento  Soil  Moisture  Accounting  (SAC-SMA)  hydrologic  model 
deployment  (Vmgt  et  al.  2009). 

Bimodal  mixed  normal  target  distribution 

Numerical  tests  to  validate  the  noted  independent  DE-MC  implementation  were  perfonned  with  a 
known  bimodal  target  composed  of  the  sum  of  two  normal  distributions,  which  is  documented  to 
be  a  difficult  problem  for  MCMC  simulation  in  that  the  distance  between  the  two  modes  makes 
jump  proposals  from  one  mode  to  the  other  problematic  (ter  Braak  2006;  Vrugt  et  al.  2009).  In 
particular,  numerical  experiments  were  perfonned  for  the  specified  known  bimodal  target: 

^(x)  =  |jVJ(-5,J„)  +  |jVJ(5>Jd)  (2) 

where  n,  x,  Nd,  5,  and  Id  represent  the  target  distribution,  a  given  chain  in  the  population,  a 
multivariate  normal  distribution  of  dimension  d,  a  vector  consisting  of  a  5  in  each  of  the  d- 
dimensions,  and  a  ^-dimensional  identity  matrix,  respectively. 

Consistent  with  ter  Braak  (2006),  all  of  the  tests  with  the  known  bimodal  target  specified  in 
Equation  2,  involved  (1)  a  predetermined  bum-in  of  1000  generations,  (2)  two  different  initial 
distributions;  viz.,  a  narrow  initial  distribution  given  by  N(0,ld)  and  a  broad  initial  distribution 
given  by  N( 2.5,25  Id),  (3)  two  values  specified  for  d  and  N;  viz.,  5  and  10,  and  100  and  1000, 
respectively,  and  (4)  y  set  equal  to  2.38/2  .  In  each  case,  100  independent  trials  were  performed 
and  uniform  random  sampling  was  employed  to  initialize  the  population.  The  subsequent 
monitoring  period  consisted  of  1000  generations.  Table  1  summarizes  the  trial  runs.  Figures  1  - 
7  provide  an  additional  summary  of  the  numerical  experiments.  For  d  =  10  and  with  the  default 
value  specified  for  y,  it  was  observed,  as  is  clearly  evident  upon  inspection  of  Figures  4  and  6, 
that  DE-MC  had  difficulty  converging  to  the  known  target  distribution.  However,  activation  of 
the  adaption  functionality  for  y,  previously  mentioned  and  briefly  discussed  above,  remedied  the 
observed  problem,  as  indicated  in  Table  1  and  also  in  Figures  5  and  7.  For  d  =  10  and  N  =  1000, 
additional  numerical  experiments,  each  consisting  of  100  independent  trials,  were  also  performed 
wherein  y  was  effectively  set  equal  to  one  every  fifth,  fourth,  and  third  generation,  in  attempts  to 
examine  how  the  update  frequency  would  impact  the  acceptance  rate  and  also  the  RMSE  for  the 
estimate  of  the  expected  value  (1.666).  Table  2  summarizes  the  results  from  these  additional 
numerical  experiments. 
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Table  1.  Summary  of  numerical  experiments  with  independent  DE-MC 
implementation  with  the  known  bimodal  target  distribution  (If  yes  for  adaption, 
then,  as  indicated  in  the  text  of  the  note,  y  is  effectively  set  equal  to  one  every 
tenth  generation;  RMSE  for  estimate  of  the  expected  value,  which  is  1.666). 

Results  based  on  the  100  independent  trials. 

d 

N 

Initial  Population 

Adaption 

Acceptance  Rate 

RMSE 

5 

100 

Narrow 

No 

0.163 

0.37 

5 

1000 

Narrow 

No 

0.157 

0.112 

5 

100 

Broad 

No 

0.166 

0.352 

5 

1000 

Broad 

No 

0.159 

0.108 

10 

1000 

Narrow 

No 

0.13 

1.609 

10 

1000 

Narrow 

Yes 

0.129 

0.124 

10 

1000 

Broad 

No 

0.199 

2.054 

10 

1000 

Broad 

Yes 

0.145 

0.23 

Figure  1.  Estimated  marginal  posterior  probability  distribution  of  x1  from  trial  33,  based 
on  using  d  =  5,  N  =  100,  and  the  narrow  initial  population.  The  solid  black 
line  depicts  the  true  bimodal  target  distribution. 
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Figure  2.  Estimated  marginal  posterior  probability  distribution  of  x1  from  trial  2,  based  on 
using  d  =  5,  N  =  1 000,  and  the  narrow  initial  population.  The  solid  black  line 
depicts  the  true  bimodal  target  distribution. 


Figure  3.  Estimated  marginal  posterior  probability  distribution  of  x3  from  trial  2,  based  on 
using  d  =  5,  N  =  1 000,  and  the  broad  initial  population.  The  solid  black  line 
depicts  the  true  bimodal  target  distribution. 
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Figure  4.  Estimated  marginal  posterior  probability  distribution  of  x5  from  trial  3,  based  on 
using  d  =  1 0,  N  =  1 000,  the  narrow  initial  population,  and  not  adapting  y.  The 
solid  black  line  depicts  the  true  bimodal  target  distribution. 


Figure  5.  Estimated  marginal  posterior  probability  distribution  of  x5  from  trial  5,  based  on 
using  d  =  1 0,  N  =  1 000,  the  narrow  initial  population,  and  setting  y  equal  to  one 
effectively  every  tenth  generation.  The  solid  black  line  depicts  the  true  bimodal 
target  distribution. 
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Figure  6.  Estimated  marginal  posterior  probability  distribution  of  x5  from  trial  4,  based  on 
using  d  =  1 0,  N  =  1 000,  the  broad  initial  population,  and  not  adapting  y.  The 
solid  black  line  depicts  the  true  bimodal  target  distribution. 


Figure  7.  Estimated  marginal  posterior  probability  distribution  of  x5  from  trial  5,  based  on 
using  d  =  1 0,  N  =  1 000,  the  broad  initial  population,  and  setting  y  equal  to  one 
effectively  every  tenth  generation.  The  solid  black  line  depicts  the  true  bimodal 
target  distribution. 
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Table  2.  Summary  of  the  impact  of  the  update  frequency  for  y. 

d 

N 

Initial  Population 

Adaption  Update 
Frequency 

Acceptance 

Rate 

RMSE 

10 

1000 

Narrow 

0.1 

0.129 

0.124 

10 

1000 

Narrow 

0.2 

0.120 

0.091 

10 

1000 

Narrow 

0.25 

0.115 

0.085 

10 

1000 

Narrow 

0.333 

0.107 

0.086 

10 

1000 

Broad 

0.1 

0.145 

0.230 

10 

1000 

Broad 

0.2 

0.127 

0.109 

10 

1000 

Broad 

0.25 

0.119 

0.093 

10 

1000 

Broad 

0.333 

0.109 

0.085  i 

SAC-SMA  hydrology  model 

Vrugt  et  al.  (2009)  described  the  DREAM  MCMC  sampler  and  compared  it  with  other  MCMC 
samplers  via  several  case  studies,  one  which  involved  inferring  the  posterior  probability 
distribution  for  thirteen  Sacramento  Soil  Moisture  Accounting  (SAC-SMA)  hydrologic  model 
parameters  for  a  SAC-SMA  deployment  to  the  1944  km“  Leaf  River  watershed  near  Collins,  MS 
using  two  years  of  mean  daily  flow  data.  The  SAC-SMA  hydrologic  model  is  used  by  the  National 
Weather  Service  (NWS)  for  flood  forecasting  throughout  the  United  States.  While  it  has  sixteen 
parameters  that  need  to  be  specified  by  the  user,  consistent  with  previous  work  (Vrugt  et  al.  2003 
and  references  cited  therein),  thirteen  were  specified  as  adjustable.  The  prior  uncertainty  ranges  of 
the  specified  adjustable  SAC-SMA  hydrologic  model  parameters  are  defined  in  Table  3. 


Table  3.  SAC-SMA  model  parameters, 
including  their  prior,  and  units. 

Parameter 

Minimum 

Maximum 

Unit 

UZTWM 

1 

150 

[mm] 

UZFWM 

1 

150 

[mm] 

UZK 

0.1 

0.5 

day-1 

PCTIM 

0 

0.1 

[-]  1 

ADIMP 

0 

0.4 

[-] 

ZPERC 

1 

250 

[-]  1 

REXP 

1 

5 

[-] 

LZTWM 

1 

500 

[mm] 

LZFSM 

1 

1000 

[mm] 

LZFPM 

1 

1000 

[mm] 

LZSK 

0.01 

0.25 

day-1 

LZPK 

0.0001 

0.025 

day-1 

PFREE 

0 

0.6 

H  1 

The  reader  is  referred  to  Vrugt  et  al.  2003,  and  references  cited  therein,  for  comprehensive 
discussions  regarding  the  SAC-SMA  hydrologic  model,  the  Leaf  River  watershed  and  its  related 
hydrologic  data  (viz.,  mean  areal  precipitation  (mm/day),  potential  evapotranspiration  (mm/day), 
and  streamflow  (m3/s))  which  was  used  to  support  the  inference  of  the  posterior  distribution  of 
the  SAC-SMA  specified  adjustable  model  parameters. 
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Consistent  with  Vrugt  et  al.  (2009),  for  the  series  of  implementation  validation  and  evaluation 
tests  employed  with  this  “real  world”  case  study  application,  squared  deviation  likelihood  was 
utilized  as  follows: 


*■(*)=  EL(s’-(*Mz 


(3) 


A 

where  S{ ,  S),  and  T  represent  the  streamflow  observations,  their  model  simulated  counterparts, 

and  the  number  of  data  comparisons  used  to  evaluate  the  likelihood  (viz.  731).  With  each  SAC- 
SMA  forward  model  call,  approximately  two  months  of  simulation  occurs  prior  to  the  window 
for  performing  comparisons  of  model  simulated  flows  with  their  observed  counterparts. 

For  purposes  of  illustration,  Figure  8  displays  results  associated  with  five  unique  realizations  of  the 
DE-MC  method  to  calibrate  the  SAC-SMA  Leaf  River  hydrologic  model.  In  particular,  it  includes 
five  plots  of  the  RMSE  computed  each  evolution.  For  each  of  the  five  DE-MC  simulation  runs, 
N=  26  and  y  was  effectively  set  equal  to  one  every  fifth  generation.  The  observed  differences  in 
Figure  8  across  the  five  RMSE  plots  is  a  function  of  the  initial  seed  specified  to  the  DE-MC 
implementation  as  all  of  the  other  input  settings  remained  constant.  The  five  plots  in  Figure  8  not 
only  underscore  the  well  known  documented  difficulty  of  calibrating  the  SAC-SMA  conceptual 
hydrologic  model  (Gupta  et  al.  2003;  Duan  et  al.  1992),  but  also  the  variability  in  the  perfonnance 
of  the  DE-MC  method.  One  can  clearly  observe  six  unique  optima  upon  inspection  of  the  five 
RMSE  plots  in  Figure  8.  Moreover,  one  can  also  observe  from  the  five  plots  that  the  number  of 
forward  model  calls  necessary  to  achieve  burn-in  for  the  DE-MC  sampler  can  exhibit  a  fair  degree 
of  variability  from  one  simulation  run  to  another.  With  the  independent  DE-MC  implementation,  a 
minimum  RMSE  value  of  13.25  was  achieved,  as  was  also  reported  upon  in  Vrugt  et  al.  (2009). 
Figure  9  is  a  trace  plot  of  the  specified  SAC-SMA  adjustable  model  parameter  UZTWM  associated 
with  the  first  chain  for  one  of  the  five  DE-MC  simulations.  Figure  10  is  a  plot  of  the  95  percent 
predictive  uncertainty  bounds  associated  with  one  of  the  five  DE-MC  simulation  runs.  The 
computed  uncertainty  bounds  are  solely  a  function  of  the  model  parameter  uncertainty  obtained  via 
sampling  from  the  posterior  probability  density  during  the  post  bum-in  monitoring  period  which 
was  arbitrarily  set  to  be  equal  in  length  to  the  number  of  simulations  required  to  achieve  DE-MC 
sampler  bum-in. 

DISCUSSION  AND  CONCLUSIONS:  This  technical  note  briefly  describes  the  Differential 
Evolution  Markov  Chain  MCMC  sampler,  which  was  identified  in  Skahill  (2012)  to  be  a  state- 
of-the-art  methodology  for  quantifying  hydrologic  model  parameter  and  predictive  uncertainty  in 
support  of  risk-based  hydrologic  design,  and  evaluation  of  project  alternatives.  This  document 
also  describes  some  of  the  salient  features  associated  with  an  independent  implementation  of  the 
DE-MC  method;  viz.,  approaches  to  initialization  of  the  population  to  be  evolved,  treatment  of 
jump  proposal  dimensions  that  are  out  of  bounds,  and  a  hybrid,  heuristic,  semi-automated 
approach  for  assessing  convergence  of  the  DE-MC  sampler.  In  addition,  this  technical  note 
summarizes  two  case  study  applications  of  the  independent  DE-MC  implementation. 

The  first  case  study  involved  a  bimodal  mixed  normal  target  distribution,  a  non-trivial 
verification  problem  in  that  in  the  past  been  demonstrated  to  be  difficult  to  resolve  for 


10 


ERDC/CHL  CHETN-VIII-7 
September  2013 


Figure  8.  A  plot  of  mean  RMSE  for  each  generation  for  five  unique  DE-MC  runs. 


Figure  9.  A  plot  of  the  evolution  of  the  first  chain  for  the  SAC-SMA  parameter  UZTWM. 
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Figure  10. 95  percent  predictive  uncertainty  envelope  which  is  solely  a  function 
of  model  parameter  uncertainty. 


classical  random  walk  metropolis  MCMC  samplers.  However,  systematic  adaption  of  one  of  the 
two  DE-MC  algorithm  input  parameters  ensured  convergence  to  the  bimodal  target  distribution. 
The  second  case  study  involved  application  of  the  independent  implementation  of  the  DE-MC 
method  to  calibrate  a  Sacramento  Soil  Moisture  Accounting  hydrology  model  deployment  using 
two  years  of  mean  daily  flow  data.  This  validation  and  evaluation  case  study  demonstrated  that  the 
DE-MC  method  is  able  to  handle  complex  response  surfaces  and  find  the  documented  minimum. 
And  that  by  sampling  from  the  equilibrium  posterior  probability  distribution  for  the  specified 
adjustable  model  parameters  one  can  then  generate  model  predictive  uncertainty  bounds  which  can 
support  risk-based  analyses.  Although,  in  the  second  case  study  example  it  was  also  illustrated  that 
the  DE-MC  method  has  the  capacity  to  exhibit  a  fair  degree  of  random  variability,  when  applied,  in 
terms  of  the  number  of  forward  model  calls  required  to  achieve  burn-in  for  the  MCMC  sampler. 

The  work  documented  herein  is  just  the  beginning  with  respect  to  research  and  development 
directed  to  the  development  of  a  state-of-the-art  robust  methodology  for  computing  model 
parameter  and  predictive  uncertainty.  Opportunities  for  future  DE-MC  related  study  and 
exploration  include: 

1.  Refine  the  current  hybrid,  heuristic,  semi-automated  approach  for  assessing  DE-MC 
convergence  by  addressing  the  question  “Can  MCMC  bum-in  be  completely  automated  and 
remain  reliable?”. 

2.  Ensure  implementation  of  the  DE-MC  algorithm  is  as  efficient  as  possible  and  reduce  the 
observed  variability  in  DE-MC  algorithm  performance. 

a.  Would  dynamic  adjustment  of  the  DE-MC  parameter  y  improve  algorithm  performance? 
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b.  Explore  applying  those  aspects  of  the  DREAM  MCMC  algorithm  (Vrugt  et  al.  2008a, 
2008b,  2009)  which  do  not  complicate  its  implementation  and  practical  usage,  such  as 
dynamic  chain  removal. 

3.  Bayesian  MCMC  provides  one  with  a  robust  means  for  assessing  model  uncertainty; 
however,  it  is  computationally  costly.  Regardless  of  the  efficiencies  derived,  port  the  final 
independent  Bayesian  MCMC  implementation  to  the  HPC  framework. 

4.  Integrate  the  independent  DE-MC  implementation  into  the  USACE-ERDC  GSSHA 
hydrology  model  (Downer  and  Ogden  2003a,  b)  as  an  alternate  simulation  mode. 

5.  Enhance  the  DE-MC  GSSHA  alternate  simulation  mode  to  also  account  for  precipitation 
forcing  uncertainty  via  specification  of  event  multipliers  similar  to  the  approach  described  in 
Vrugt  et  al.  (2008b). 

6.  Modify  the  independent  DE-MC  implementation  to  account  for  model  structural  uncertainty 
via  specification  of  a  likelihood  measure  that  accommodates  for  the  treatment  of  serial 
correlation  of  the  residuals.  Explore  alternative  means  to  quantify  model  structural 
uncertainty. 

7.  Incorporate  into  the  independent  DE-MC  implementation  automated  means  for  computing 
the  necessary  length  of  the  post  burn-in  monitoring  period. 

8.  Provide  a  more  complete  treatment  and  discussion  of  the  various  options/features  of  the 
independent  implementation  of  the  DE-MC  method  to  provide  potential  users  with  a  clear 
picture  regarding  recommended  usage  of  various  options  such  as  use  of  URS  versus  LHS 
for  initializing  the  population  to  be  evolved  or  the  treatment  of  jump  proposal  dimensions 
that  are  out  of  bounds. 

9.  Clearly  demonstrate  the  capacity  of  the  Gelman  and  Rubin  diagnostic  to  prematurely  assess 
MCMC  sampler  convergence. 

10.  Perfonn  comprehensive  case  studies  that  demonstrate  use  of  the  state-of-the-art  model 
uncertainty  analysis  method  to  support  practical  risk-based  hydrologic  and  environmental 
application  settings. 

ADDITIONAL  INFORMATION:  For  additional  information,  contact  Dr.  Brian  E.  Skahill  at 

Brian. E. Skahill(a),usace. arm y. mil.  This  CHETN  should  be  cited  as  follows: 

Skahill,  B.  E.  (2013).  Initial  development  and  testing  of  a  state-of-the-art  method 
to  quantify  hydrologic  model  uncertainty.  ERDC/CHL  CHETN-VIII-7. 
Vicksburg,  MS:  US  Army  Engineer  Research  and  Development  Center.  An 
electronic  copy  of  this  CHETN  is  available  from  http: //chi  erdc.  usace.army. 
mil/chetn. 
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